Rotating star initial data for a constrained scheme 
in numerical relativity 



O 
O 

(N 



cr 



Lap-Ming Lin and Jerome Novak 

Laboratoire de l'Univers et de ses Theories, Observatoire de Paris, F-92195 
Meudon Cedex, France 



Abstract. A new numerical code for computing stationary axisymmetric 
rapidly rotating stars in general relativity is presented. The formulation is based 
on a fully constrained-evolution scheme for 3 + 1 numerical relativity using the 
Dirac gauge and maximal slicing. We use both the polytropic and MIT bag model 
equations of state to demonstrate that the code can construct rapidly rotating 
C ~ ) ' neutron star and strange star models. We compare numerical models obtained by 

ff) ' our code and a well-established code, which uses a different gauge condition, and 

show that the two codes agree to high accuracy. 
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In the 3 + 1 formalism of general relativity, the Einstein equations are decomposed 
into a set of four constraint equations and six evolution equations ^12- Solving the 
(elliptic-type) constraint equations at each time-step in multidimensional simulations 
is in general not feasible as it is computationally expensive for most numerical 
techniques. Hence, a free-evolution approach, i.e., solving the constraint equations 
only for the initial data and performing the evolutions without enforcing the 
constraints, is generally favoured over a fully constrained scheme (solving the 
constraints at each time step) in three-dimensional numerical simulations. While 
mathematically the constraints are preserved by the evolution equations, in practice 
small constraint violations due to numerical errors typically grow quickly to a 
significant level that make the solution unphysical and plague the simulations. 
Although numerous techniques to control the growth of constraint violations have 
been developed (e.g., EJ 03 El 13 E])j it is not clear yet to what extent they can 
control the constraint violations successfully in general. 

Fully constrained-evolution scheme has been used in the past only in spherically 
symmetric or axisymmetric problems (e.g., 0E31E1E])- The main advantage of a 
fully constrained scheme is that the constraints are fulfilled to within the discretisation 
errors, and the constraint-violating modes do not exist by construction. Recently, a 
new formulation for 3 + 1 numerical relativity based on a fully constrained-evolution 
scheme is proposed by Bonazzola et al. |13| . The chosen coordinate conditions 
(maximal slicing and Dirac gauge, see Sec. 12.3(1 . together with the use of spherical 
components of tensor fields, reduce the ten Einstein equations to a system of five quasi- 
linear elliptic equations, which are solved by efficient multi-domain spectral methods, 
and two quasi-linear scalar wave equations |13|. The Dirac gauge is used to fix the 
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remaining three degrees of freedom. The stability of the proposed scheme has been 
demonstrated for a three-dimensional pure gravitational wave spacetime ^3] ■ 

In the proposed constrained scheme, the coordinates are fixed by the Dirac gauge 
on each hypersurface (with maximal slicing condition), including the initial one. This 
implies that initial data must be prepared in the same coordinate choices in order 
to perform dynamical evolutions. As an advantage of this gauge- fixing, stationary 
solutions of the Einstein equations can be computed within the same framework 
simply by setting the time derivative terms to zero in the equations. The aim of 
the present work is to construct rotating-star initial data for this new formulation of 
3 + 1 numerical relativity. For this purpose, we have developed a numerical code to 
calculate stationary axisymmctric rotating star models based on the Dirac gauge and 
maximal slicing. 

The purpose of this paper is to present the formulation of the problem and 
the tests that have been done to validate our numerical code. The numerical code 
can be used to provide initial rotating star models for hydrodynamics simulations 
in full general relativity within the new formulation. Our emphasis here is put 
on the comparison of the accuracy between our code and a well-established code 
LORENE-rotstar [141 1151 H?)| . In particular, we demonstrate that our numerical code 
can compute rapidly rotating neutron star and strange star models to high accuracy. 
Unless otherwise noted, we use units such that G = c = 1. Latin (Greek) indices go 
from 1 to 3 (0 to 3). 

2. Formulation 

2.1. The 3 + 1 decomposition 

In this section, we give a brief description of the 3+1 formulation of the Einstein 
equations in order to define our notations (see, e.g., P3E] for details). In the 3 + 1 
formalism, the spacetime is foliated by a family of spacelike hypersurfaces S t , labelled 
by the time coordinate t. Introducing a coordinate system (x l ) on each hypersurface, 
the line element may be written as 

ds 2 = -N 2 dt 2 + 7y(dx* + P l dt){dx 3 + ftdt), (1) 

where N is the lapse function, (3 l is the shift vector, and 7^ is the 3-metric induced 
by the spacetime metric g a p onto each hypersurface E t 

7q/3 : = 9aj3 + n a np. (2) 

Here n a := —NV a t is the unit normal to S t , where V Q is the covariant derivative 
associated with the spacetime metric g a p. The stress-energy tensor T a/3 is decomposed 
as 

T a f } = En a n l3 + n a J + J a n + S a(} , (3) 

where E := T a pn a n^, J a := — ^J^T^n" , and S a p :— 7 Q f7 | g / T' AH , are the energy density, 
momentum density, and the stress tensor, as measured by the observers of 4-velocity 
n a (the so-called Eulerian observers O e ). 

The evolution of the 3-metric 7^ is governed by 

^7«-^/97« = -2^ 1 (4) 
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where C is the Lie derivative operator and K t j is the extrinsic curvature of E t . The 
evolution equation for is 

Tjy A., - CpKn = - DiDjN + .V {/,',, - , A '•*' + A'A, ; 

+ 47r[(5-^) 7ij --2^]}, (5) 

where is the covariant derivative associated with the 3-metric 7y , i?^ is the Ricci 
tensor associated with this 3-metric, K := A' l j is the trace of the extrinsic curvature, 
and S :— S 1 ^ In the 3+1 formulation, the full set of Einstein equations is equivalent 
to the above evolution equations, together with the Hamiltonian constraint 

R + K 2 - K l3 K ij = 16ttE, (6) 
and the momentum constraint 

DjKf - DiK = 8-jrJi, (7) 
where R := R l i is the three-dimensional Ricci scalar. 



2.2. The matter sources: uniformly rotating fluid 

We assume that the matter consists of a perfect fluid with a stress-energy tensor 

T a/3 = (e + P)u a vP + Pg a/3 , (8) 

where u a is the 4- velocity of the fluid; e and P are respectively the energy density and 
pressure, as measured by the the fluid comoving observer Of. Defining the Lorentz 
factor linking the two observers O e and Of by 

T := ~n a u a = Nu\ (9) 
the energy density E in Eq. © can be written as 

E = T 2 (e + P) - P. (10) 
The momentum density J 1 is 

,p = (e + p)v\ (ii) 

where the fluid 3- velocity v l is related to the spatial components of the fluid 4- velocity 
v? by 



N 



(12) 



Note that the Lorentz factor can be expressed as L = (1 — u 2 ) -1 / 2 , where v :— (viV 1 ) 1 / 2 
is the "physical" fluid speed as measured by the Eulerian observer O e . The stress 
tensor Sij is given by 

S ij = (E + P)v i v j +P'Yi j . (13) 

Up to this point, we have not made any assumptions on the spacetime and 
fluid flow. Now we consider that the spacetime is stationary, axisymmetric, and 
asymptotically flat. These assumptions imply the existence of two Killing vector fields: 
£, which is timelike at spatial infinity; x, which is spacelike everywhere and with closed 
orbits. Furthermore, the Killing vectors commute ^7], and hence we can choose an 
adapted coordinate system (t, x 1 , x 2 , ip), such that £ = d/dt and x = d/dip are the 
coordinate vector fields (see also J3]). We choose the remaining two coordinates to 
be of spherical type (i.e., x 1 — r and x 2 = 8). 



Rotating star initial data for a constrained scheme in numerical relativity 



4 



We also impose the so-called circularity condition on the stress-energy tensor: 

= K a + "X a , (14) 

r%-x? = \c + <rx a - (is) 

This condition is equivalent to the absence of meridional convective currents, and 
implies that the fluid 4-velocity is given by 

*=«*f4 + n £V (is) 



^dt dip 

where f2 := u v /u* is the fluid coordinate angular velocity, and can be interpreted as 
the fluid angular velocity as seen by an inertial observer at rest at infinity. A theorem 
of Carter JSj shows that the circularity condition implies that the line element can 
be written as 

ds 2 = -N 2 dt 2 + j vv (dtp + f3 v dt) 2 + lrr dr 2 + 2 lr9 drd0 + ^edO 2 .{17) 

Notice that only the (^-component of the shift vector is nonzero and we have not 
specified the gauge choice at this point. The so-called quasi-isotropic gauge (see 
Sec. corresponds to j r e — 0, fee = r 2 j rr , while the Dirac gauge relates the 
metric components by differential equations (see Eq. H26[) ~). 

The equation of stationary motion follows from the projection of the conservation 
equation V Q T Q ^ = normal to the 4-velocity u a . In this paper, we focus on the 
case of uniformly rotating star (i.e., Q, is a constant). We also assume that the 
fluid is barotropic. In this case, the equation of stationary motion can be integrated 
analytically and is given by (see, e.g., |14p 

iJ + lniV-lnr = const., (18) 

where H is the log-enthalpy defined by 

2.3. The metric equations 

Here we summarise the full set of Einstein equations in the constrained-evolution 
scheme based on the Dirac gauge and maximal slicing. We refer the reader to Ref. 13 
for details and derivations. First, we define a conformal metric jij by 

7„ (20) 

with the conformal factor ^> defined by 

where fij is a flat 3-metric, given by the asymptotic condition on 7y . The four 
constraint equations © and (JJJ, together with the maximal slicing condition K = 0, 
result in two scalar equations for the lapse and conformal factor; and one vectorial 
elliptic equations for the shift vector. 
The lapse function N is given by 



AN = ^ 4 N 



4tt{E + S) + A kl A kl - h kl V k V[N - 2D k $D k N, (22) 



,kl 
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where T>i is the covariant derivative associated with the flat metric /y and its 
contravariant component is defined by T> % :— f i3 T>f, A := f %3 T> i 'Dj is the flat-space 
Laplacian operator; Di is the covariant derivative associated with the conformal metric 
jij and its contravariant component is D l := 7 y Dj (with the inverse conformal metric 

part of the conformal extrinsic curvature A 13 is defined by 

1 



defined by Jik7 = <V)- The quantity is defined by $ :— In "J". The traceless 



* 4 K %3 



-f 3 K 



while the tensor field Aij is defined by 



A; 



likljiA 



K, 



(23) 



(24) 



The tensor field h 13 on the right-hand side (RHS) of Eq. (|22|l is the deviation of the 
inverse conformal metric from the inverse flat metric defined by 

h i3.-f3_fj. (25) 

In the proposed constrained scheme, the Dirac gauge condition is given by 

Vjh i3 ' = 0. (26) 

Next, the conformal factor '5 (or equivalently Q := ^ 2 N) is determined from 

3 7 
4^ 



AQ = ~ h kl V k ViQ + V 6 



N 4ttS + 



:A kl A kl 



2^ 



N ( T + ® h ®® k ® 



D k <5>D k N 



where the quantity i?* on the RHS is given by 

R* := i7 fc ^fe^ m "A7mn - \l U V k h mn V n ^ ml . 

The elliptic equation for the shift vector is 
f 



(27) 



(28) 



A/3 1 



-V\Vjf3 J ) = 16tt7V* 4 J 4 + 2A y X>jiV - 12iVA lJ D J $ - 2NA l kl A 



ki 



I 



-h kl V k V^ - -h lK V k V0 



where the tensor field A". . is defined by 



A* 



l u {Vijij+Vjja-Vijij). 



(29) 



(30) 



Now we turn to the dynamical part of the Einstein equations. In the proposed 
constrained scheme, one solves for the tensor field h lJ instead of ^f 3 . The evolution 
equation in this formulation is given by a flat-space (second-order) wave equation for 
h %3 (see Eq. (85) of Ref. pi])- As mentioned in Seed one advantage of using the 
Dirac gauge to fix the coordinates on each slice S t is that stationary solutions of 
the Einstein equations can be computed within the same scheme, simply by setting 
the time-derivative terms to zeros in the equations. This is possible because of the 
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existence of the Killing vector £ = d/dt as discussed in Sec. 12.21 This reduces the 
wave equation for h lJ to the following tensorial Poisson-like equation: 

Ah ij = ^ IcpCptii + ~V k p k £ph ij + ^V k Q (V l h jk + V j h lk - V k h lj ) - 2NS lj 



£(3(V k/ 3 k ) + -{V k p k f 



h IJ 



(31) 

where the notation {Lffy 3 stands for the conformal Killing operator associated with 
the flat metric acting on the shift vector 

(L0yi := V l /3 3 + V 3 p l - 2 -V k p k p. (32) 

The tensor field S l] on the RHS of Eq. (gTJ is given by 



S 13 = 



i r 

3 



N 



R*+8D k <I>D k <I> 



8D k <PD k Nf- 



+2N 



i 



yuA^A" - 4tt ( * 4 5 y - -Sf 3 



(33) 



with 
= 



j k n mn v m h lk v n h 3 ' 



+lnlV k h mn (f k V m h 31 + y 3k V m h 



-7 7 J D k h Di^ n 



(34) 



Furthermore, after setting the time-derivative term of h l ° to zero, the conformal 
extrinsic curvature A 13 as defined in Eq. (|23(l is deduced from (see Eq. (92) of |13p 



A i3 = — 
2N 



(Lf3y 3 - Cph? 3 - -V k f3 k W 



(35) 



In summary, to calculate stationary axisymmetric uniformly rotating star models 
in the framework of the constrained scheme ^H] based on the Dirac gauge and maximal 
slicing, one needs to solve for two scalar elliptic equations <|22[) and l|27|) respectively for 
N and 'J, a vectorial Poisson-like equation (|29|l for and a tensorial elliptic equation 
13111 for h lJ , together with the first integral of motion Eq. (|18|) for the matter. The 
numerical procedure on how to solve this system of equations is described in Sec. 13.11 



2.4- Global quantities 

We list here various global quantities relevant to axisymmetric rotating-star 
spacetimes. These gauge invariant quantities are useful to estimate the accuracy 
of our numerical code as they provide a direct comparison between our code and a 
different code, which uses a different gauge condition, as presented in Sec. 13.31 
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Given a baryon current nu a , where n is the number density in the fluid frame, 
the baryon mass of the star is expressed as 

Mb = m B J [-n a (nu a )} dV = m B / TndV, (36) 

where me is the baryon rest mass, dV — ^fd 3 x is the proper 3-volume element (with 
7 being the determinant of the 3- metric), and we have used Eq. @ in the second 
equality. As we follow Ref. J3| to expand all tensor fields onto the spherical basis 
( e i) = iff i r s ln e ~§ip ) ' w hich is orthonormal with respect to the flat metric, the 
proper volume clement is written explicitly as dV = a/tV 2 sin 8 'drdO 'dip. Notice that 
we denote by j the determinant of the 3-metric expanded onto the basis (ej). Here 
and afterward we denote by (f, 9, ip) indices of specific components on the orthonormal 
basis (ej). 

The gravitational mass M g is given by the Komar integral (see Eq. (11.2.10) of 

PI) 

M g = 2 J (r afj - \T\g al \ n a C dV : (37) 
where ( a is the timelike Killing vector discussed in Sec. 12.21 Explicitly, we have 

M g = J [N(E + S) - 2J^\ dV, (38) 
The total angular momentum J is given by (see Problem 6 of jiy] ') 

J = - I T al3 n a X ' 3 dV= [ J^rsinedV, (39) 



where \ a is the axial Killing vector of the spacetime. The rotational kinetic energy 
for a uniformly rotating star is T = \QJ ■ The gravitational potential energy is 

W = M P + T-M g , (40) 

where M p is the proper mass of the star defined by 

M p = J [-n a (eu a )] dV = J TedV. (41) 

Note that W is defined to be positive. 

Furthermore, two relativistic virial identities (the so-called GRV2 and GRV3) 
have been proved to be useful for checking the consistency and accuracy of numerical 
solutions of rotating relativistic star models. The three-dimensional virial identity 
GRV3 20] is a relativistic generalisation of the Newtonian virial identity, valid for any 
stationary and asymptotically flat spacetime. The two-dimensional virial identity 
GRV2 |23 is valid for any asymptotically flat spacetime (without any symmetry 
assumption). The two virial identities are integral relations between the matter 
and metric fields (see [20112] for expressions). In practice, we define the quantities 
GRV2 := 1 1 — Aa J and GRV'i := 1 1 — A3 1 as the error indicators for the virial identities, 
where A2 and A3 are defined via the integral relations, such that exact solutions of the 
Einstein equations satisfy GRV2 = GRV3 = (see [221 )■ Note that these identities 
are not imposed during the numerical calculation, and hence are useful indicators for 
checking the accuracy of numerical results. 
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3. Numerical results 

3.1. Numerical procedure 

To calculate stationary axisymmetric rotating stars within the Dirac gauge and 
maximal slicing, we solve the nonlinear elliptic equations described in Sec. 12.31 
iteratively by means of multi-domain spectral methods |2S1 123 m spherical 
coordinates. The code is constructed upon the LORENE C++ library We use 
three spherical numerical domains to cover the whole hypersurface E t . Specifically, 
we use one domain to cover the whole star and one domain for the space around the 
star (typically to about twice the stellar radius). The last domain covers the space 
out to spatial infinity by means of a compactification u = 1/r |14| . 

In each domain, we use Ng collocation points in the polar direction and N v = 1 
point in the azimuthal direction for the spectral method. For the radial direction, we 
can choose to have different numbers of collocation points in different domains. In 
Sec. 13.31 we use the notation N r = (N r %, N r 2, N r a), where N r i denotes the number of 
points in the first domain etc, to specify the grid structure in the radial direction. 

The numerical iteration procedure is briefly described here. For a given equation 
of state, we choose fi and the central value of the log-enthalpy H (see Eq. tE^l ) as the 
physical parameters that specify the rotating star model. First we start with an initial 
guess by setting all the metric quantities to their flat spacetime values, together with 
a spherically symmetric distribution for the matter sources. The iteration procedure 
begins by solving Eqs. (122[l and (|29|l respectively for the corresponding lapse and shift. 
Thus, we obtain the only nonzero component of the 3- velocity = {Sir sin 9 + j3 v ) /N 
(see Eq. (|12|) V and hence the Lorentz factor T. Next we use the first integral of motion 
(Eq. I|18jl) to obtain H, from which we deduce the pressure P and the energy density 
e through the EOS. Finally, we solve Eqs. l(77|) and (|3"T)l respectively for Q and h l K 
The iteration procedure continues until the relative difference in H throughout the 
whole star between two consecutive steps is smaller than some prescribed value. 

The resolutions of the scalar Poisson equations for N and Q, and the vectorial 
elliptic equation for /3* have been described in details in • The technique for solving 
the tensorial Poisson equation (|31|) is described in | Appendix A| 



3. 2. Equation of state 

In Sec. 13.31 we present various tests that have been done to validate our numerical 
code. For this purpose, we use a polytropic EOS in the following form to construct 
rotating neutron star models: 

P = «n 7 , (42) 

where n and 7 are constants. The number density n is related to the energy density 
e by 

e = m^n H n 7 , (43) 

7- 1 

where the baryon mass me = 1.66 x 10~ 24 g. In particular, we take 7 = 2 and 
k — 0.03p nuc c 2 /n 2 uc , where p nuc = 1.66 x 10 14 g cm~ 3 and n nuc = 0.1 fm~ 3 . For this 
EOS, the log-enthalpy is given analytically by 



H = In 



1 



«7 



,7-1 



"i B (7 _ 1 ) 



(44) 



Rotating star initial data for a constrained scheme in numerical relativity 9 




20 40 60 

Iterations 



Figure 1. Convergence towards zero of the relative difference in H throughout 
the star between two consecutive steps for a non-rotating polytropic star model 
(see text). 



We also use the simplest MIT bag model EOS, with noninteracting massless 
quarks, to construct rapidly rotating strange stars. The EOS is given in the following 
form (see, e.g., 

P = \an A ' 3 - B, 

e = an 4/3 + 5, (45) 

where B is the MIT bag constant and the parameter a = 9ir 2 / 3 hc/4: = 
952.371 MeV fm. We choose B = 60 MeV fm" 3 in this work. The stellar surface 
is characterised by the properties of strange matter at zero pressure: the number 
density n = 0.28665 fur 3 and the mass density p Q = 4.2785 x 10 14 g cm -3 . The 
log-enthalpy is related to n simply by H = ln^ra/rao) 1 / 3 . The MIT bag model EOS is 
useful to test our numerical code in the highly relativistic regime, since strange stars 
can reach higher compactness ratios and rotation rates than ordinary neutron stars. 

3. 3. Tests of the numerical code 

To test our numerical code, we start with a non-rotating star modelled by the 
polytropic EOS described in Sec. 13.21 The central value of the log-enthalpy is 
H = 0.2308 (or equivalently the energy density eo = 4.889 j o nuc c 2 ) and the baryon 
mass of the star is Mb = 1.6M©. The star has a gravitational mass M g = 1.4866M Q 
and a compactness ratio M g /R = 0.147, where R is the circumferential radius. In the 
numerical calculations, we use a parameter en (typically set to be 10" 10 or smaller) to 
control the iteration procedure and the precision of the numerical models: the iteration 
(see Sec. I3.1|) is stopped if the relative difference in H throughout the whole star 
between two consecutive steps is smaller than en- In Fig. [I] we show the convergence 
of the relative difference in H towards zero with the number of iterations using radial 
collocation points N r = (33,33, 17). We see that a precision of 10~ 15 is achieved for 
the numerical result within 40 iteration steps. After that, the accuracy is limited by 
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Figure 2. Convergence towards zero of GRV2 and GRV3 with the number of 
collocation points N r \ inside the star for the same model as shown in Fig. 



the round-off errors. The solution also satisfies the virial identities to the level of 

io- 15 . 

Next, we test the convergence property of the numerical code with respect to 
increasing number of radial collocation points using the same model. In particular, we 
vary the number of points in the first numerical domain (i.e., inside the star), while 
keeping the points in the other two domains fixed to N r 2 = 33 and N r s = 17, and 
we choose e H = 10~ 15 in this test. In Fig. H we plot \og{GRV2) and log(Gi?V3) 
together against the number of points N r \ in the first domain. It is seen clearly 
that both GRV2 and GRV3 converge exponentially towards zero with the number of 
points, as expected for spectral methods. The accuracy is limited by the round-off 
errors for N r % > 30. 

Starting from the above non-rotating polytropic model, we then construct a 
sequence of increasing uniformly rotating configurations at fixed baryon mass Mb = 
1.6Mq, up to the mass-shedding limit. The accuracy of the numerical models are 
estimated by comparing the results to those obtained by a separated code LORENE- 
rotstar [141 115) . which uses the so-called quasi-isotropic gauge to construct rotating 
relativistic stars. LORENE-rotstar is a well-established code which has been tested 
extensively and compared with a few different numerical codes |22| . A comparison 
of some of the gauge invariant quantities for the sequence between our code and 
LORENE-rotstar is given in Tabled The parameters used to obtain the numerical 
results are N r — (33, 17, 17), Ng = 17, and e# = 10~ 10 for both codes. In the table, 
for each value of the rotation frequency /, we display the results obtained from our 
numerical code based on the Dirac gauge in the first row. Below this are the results 
obtained from LORENE-rotstar. Table shows that the two sets of results agree 
to high accuracy. In particular, the overall discrepancy between the two different 
codes is consistent with the errors in the virial identities, which increase with the 
rotation frequency. Note that the errors in the virial identities for the non-rotating 
configuration listed in the table is about 10~ 9 instead of 10 -15 as shown in Fig. [5] 
This is due to our choice of using e# = 10~ 10 in this test. 

We note, however, that the numerical error no longer decreases exponentially with 
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Table 1. Comparison between our numerical code based on the Dirac gauge 
(first row for each given frequency) and a well-established code LORENE-rotstar 
(second row), which uses a quasi-isotropic gauge for the coordinates, for a 
sequence of 7 = 2 polytropic neutron star models with fixed baryon mass 
Mf, = 1.6Mq. Listed are the rotation frequency /, gravitational mass M g , total 
angular momentum parameter J/Mg, ratio of the rotational kinetic energy to the 
potential energy T/W, equatorial circumferential radius R eq , and errors indicators 
in the virial identities GRV2 and GRV3. 



f (Hz) 


Mo (Mm) 


" / q 


T/W 




GRV2 


GRV3 





1.486610961 








14.91222928 


7e-9 


2e-9 




1.486610965 








14.91222929 


9e-9 


le-9 


100 


1.486837016 


0.066036124 


0.00107222950 


14.9609161 


3e-9 


3e-8 




1.486837013 


0.066036126 


0.00107222954 


14.9609165 


3e-8 


5e-8 


200 


1.487539907 


0.13421663 


0.004395974 


15.113935 


2e-7 


5e-7 




1.487539902 


0.13421665 


0.004395975 


15.113937 


3c-7 


7e-7 


300 


1.4888008 


0.2072617 


0.010340091 


15.395905 


3e-7 


6e-7 




1.4888007 


0.2072618 


0.010340095 


15.395908 


2e-8 


6e-7 


400 


1.49080352 


0.289551 


0.01973704 


15.865576 


lc-6 


3e-6 




1.49080359 


0.289550 


0.01973701 


15.865574 


7e-7 


4e-6 


500 


1.493991 


0.390430 


0.0345952 


16.67896 


4e-6 


8e-6 




1.493990 


0.390432 


0.0345954 


16.67900 


6e-6 


le-5 


550 


1.4964095 


0.455303 


0.0457202 


17.35894 


2e-5 


2e-5 




1.4964092 


0.455307 


0.0457206 


17.35898 


le-5 


2e-5 


600 


1.500054 


0.54397 


0.062368 


18.5382 


5e-6 


2e-6 




1.500055 


0.54398 


0.062369 


18.5383 


3e-6 


3e-6 


640 « h 


1.506928 


0.695855 


0.0929183 


22.1467 


2e-5 


5e-5 




1.506929 


0.695857 


0.0929188 


22.1469 


le-5 


6e-5 



the number of grid points as in the non-rotating case, but as a power-law (see Fig.^J, 
due to the discontinuities in the derivative of the matter fields at the stellar surface. 
The non-rotating model is free from any such phenomenon because the stellar surface 
is at the boundary between two spherical numerical domains. For rotating models, 
because of the flattening of the stars, the stellar surface no longer coincides with 
the boundary of the domains, and hence the spectral method loses its exponential- 
convergence property. Such difficulty associated with rotating stars can be handled by 
the adaptation of the numerical domains to the stellar surface as developed in [23 ■ We 
plan to improve our numerical code by implementing this surface-adaptation technique 
in the near futurej. Nevertheless, even without a surface-adaptation technique, Tabled 
shows that both numerical codes still agree to high accuracy and achieve a precision of 
10~ 5 for a configuration rotating near the mass-shedding limit (i.e., the / = 640 « fk 
Hz configuration in Table ^ where fk is the Kepler frequency) . To visualise the 
gravitational field generated by a rotating star, we plot in Figs. 14171 the non- vanishing- 
components for the metric field h % i (namely, h rr , h r9 , h" , and h^) for the rotating 
star model with / = 640 Hz. In these figures, we show the iso-contours of the fields in 
the meridional plane, where solid (dashed) lines indicate positive (negative) values of 
the fields. The thick solid lines represent the stellar surface. Finally, the dot-dashed 
circles represent the boundary between the first two spherical numerical domains. 
The figures show clearly that the gravitational field is dominated by the quadrupole 
moment. 

To further calibrate our numerical code against LORENE-rotstar, we now 

I Such numerical technique is already available in LORENE-rotstar as described in 1151 . However, 
in order to compare to our numerical code, the results obtained from LORENE-rotstar as listed in 
Tableware based on fixed spherical numerical domains. 
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Figure 3. Convergence behaviours of GRV2 and GRV3 for the / = 400 Hz 
rotating star model listed in Table Note that the plot is in log-log scale. The 
best-fit to the data points shows that GRV2 (GRV3) decreases as N~-^ 3 (N~^ 9 ). 



Table 2. Comparison between our numerical code (LORENE-rotstar_dirac) and 
LORENE-rotstar for a rapidly rotating strange star model. The star has a baryon 
mass Mj, = 2.2Mq, gravitational mass M g = 1.719Jl4g, and rotation frequency 



/ = 1000 Hz. 




rotstar_dirac 


rotstar 


rel. diff. 


M 9 (M Q ) 


1.7194 


1.7198 


2e-4 


J/Ml 


0.5940 


0.5945 


8e-4 


T/W 


0.0888 


0.0890 


2e-3 


R eq (km) 


12.425 


12.433 


7e-4 


GRV2 


7e-4 


2e-4 




GRV3 


le-3 


6e-4 





compare the two codes using a very relativistic and rapidly rotating strange star 
model. The matter is described by the MIT bag model as described in Sec. 13.21 This 
configuration, shown in Fig. El has a baryon mass M& = 2.2Mq and compactness ratio 
Mg/Req = 0.204 (with the gravitational mass M g = 1.719M0 and the circumferential 
equatorial radius R eq — 12.425 km). The rotation frequency is / = 1000 Hz. 
In Table [3 we show the values of various physical quantities obtained from both 
numerical codes, together with the relative difference between them. As in the case of 
the polytropic EOS model, the discrepancy of the two numerical codes is consistent 
with the errors in GRV2 and GRV3. Note also that, even with a strong density 
discontinuity at the strange star surface, our numerical model still achieves a precision 
of 10" 3 . 



4. Conclusion 



In this paper we have developed a computer code (LORENE-rotstar_dirac) to 
construct relativistic rotating stars within the framework of a new constrained- 
evolution formulation of the 3+1 Einstein equations based on the Dirac gauge and 
maximal slicing U^J. As the Dirac gauge fixes the spatial coordinates on each time 
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Figure 4. Iso-contours of the metric component h rr in the meridional plane for 
the / = 640 Hz rotating polytropic model given in Tabled The solid (dashed) 
lines indicate positive (negative) values of the field. The thick solid line represents 
the stellar surface. The dot-dashed circle is the boundary between the first two 
spherical numerical domains. 
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Figure 7. Same as Fig. 171 but for the component K** . 
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Figure 8. Meridional plane cross section of a rapidly rotating strange star. The 
star has a baryon mass Mj, = 2.3Mq, gravitational mass M g = 1.787Mq, and 
rotation frequency / = 1000 Hz. The lines are iso-contours of the log-enthalpy H. 
The thick solid line represents the stellar surface. Outside the star, H is defined 
by the first integral of motion Eq. 1181 . 

slices, including the initial one, this formulation can be used to compute stationary 
solutions of the Einstein equations simply by setting the time derivative terms in 
various equations to zeros. The system reduces to two scalar elliptic equations for 
the lapse function N and conformal factor \1/ (equivalently for Q := ^ 2 N), a vectorial 
elliptic equation for the shift vector f3 l , and a tensorial elliptic equation for a tensor 
field h % i , We couple this system of equations to the first integral of motion for the 
matter, and solve the equations iteratively using multi-domain spectral method. 

We have demonstrated that this formulation can be used to compute stationary 
rotating equilibrium configurations to high accuracy. In particular, we used the 
polytropic EOS and MIT bag model to calculate rotating neutron star and strange 
star models respectively. We compared our code to a well-established code LORENE- 
rotstar, which uses a quasi-isotropic gauge to fix the coordinates, and found that the 
global quantities of the numerical models obtained from the two codes agree to high 
accuracy. The discrepancy between the two codes is consistent to the errors in the 
virial identities. 

Finally, we remark that the proposed constrained-evolution scheme |13] is 
particular well suited to the conformally-flat relativistic hydrodynamics code, with 
a metric solver based on spectral methods and spherical coordinates, developed by 
Dimmelmeier et al. in the so-called Mariage des Maillages (MDM) project 120]. The 
numerical code that we described in this paper can be used to generate rotating-star 
initial data for hydrodynamics simulations in full general relativity within the new 
constrained-evolution scheme f° r the MDM project. 
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Appendix A. Resolution of the Poisson equations for h 1 ^ 

Here we describe the numerical strategy used to solve the tensorial Poisson equation 
13111 , imposing that the solution h 1 ^ satisfies the gauge condition l|2tj[l and be such that 
the conformal metric has a unitary determinant: 

del /'-' • h") =1. (A.l) 

Note that this relation follows directly from the definition of the conformal factor in 
the proposed constrained scheme (sec Eqs. (I2U|) and J2U), together with the condition 
det/j-- = 1 in the orthonormal basis (ej) (see Sec. 12.40 . In Ref. one would 

solve two (scalar) Poisson equations: for h rr and the potential fj, (see Eq. l|A.8|l l: the 
other four components are deduced from the three gauge conditions and the non- 
linear relation 1A.1|) through an iteration. The drawback of this method is that some 
components of h 1 ^ are calculated as second radial derivatives of h rr and \i. Since 
the source of Eq. I|31(l contains second-order radial derivatives of h l \ one needs to 
calculate fourth-order radial derivatives of h rr and fi, which are solutions of scalar- 
like Poisson equations with matter terms on the RHS. In the case of neutron stars, it 
is quite often that radial density profiles have a discontinuous derivative at the surface 
of the star. Therefore, h rr and fi admit discontinuous third-order radial derivatives 
and their fourth-order derivatives cannot be represented at all by means of spectral 
methods. A solution could be to use adaptive mapping: the boundary between two 
spectral domains coincides with the (non-spherical) surface of the star (see Still, 
the evaluation of a fourth-order radial derivative introduces too much numerical noise, 
even using spectral methods. We have therefore chosen to use a different approach, 
detailed hereafter. 

Instead of using directly all the spherical components of the tensor h tJ , we use 
only the ff- component, the trace h = fijh? 3 and the four potentials n, fi, W and X 
defined as follow, in the orthonormal basis (ej): 

r\d9 sm9 dip)' { ' 

dn dfi 
~thp~ + ~d9 
and 

_ d 2 W l_dW 1 d 2 W 9 / 1 dX \ 

d6 2 tan 6 d6 sin 2 6 «V 80 ) ' ( 

h ^PX_ _^dX _^^X d_ (J_dW\ . 

d6 2 tanfl 89 sin 2 9 dip 89 \sm9 dtp)' 1 ' ' 

with P = (h — h^)/2. These equations can be inverted to compute the potentials 
in terms of the angular part of the Laplace operator 

d 2 1 d _L £P_ 

dip ~W + i^T9d9 + ^9W ( ' 



r ^:i_L^ + |i (a.3) 
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giving 

/ dhr* h f $ 1 dhf S \ 

. , A . ttt <9 2 p 3 OP i a 2 p „„ 
A ev (A ev + 2)^= — + ——-—— -2P 

(A.9) 





) 

- + 


h fe 


1 dh^ 


^ de 




7, + 

tanfe* 


sin 9 dp 


( dh^ 


h f<P 


1 dhf S 


\ 89 




tan# 


sin 6 dp 


d 2 P 


3 


dP 


1 d 2 P 


dO 2 + 


tan 


9 89 


sin 2 9 dip 2 


2 


d 


( dh s * 




smO 


dip 


{ 09 


+ tan 6 J ' 



A ev A e „ + 2)X = — — + — — — , - 2/i^ 

cw 2 tan 90 sin 9 dip 2 

2 9 /aP P 



sinfldyj V 961 ' tan^ ' 
These quantities r),fj,,W,X are interesting for, at least, two reasons: first they 
can be expanded onto a basis of scalar spherical harmonics Y™(9,ip), which are 
often used in the framework of spectral methods, for they are eigenfunctions of the 
angular Laplace operator i|A.6|) . Furthermore, the three gauge conditions i|26|) can be 
reformulated in terms of these potentials: 

dh ff 3h ff * h 

-5— + + ^A 0lp r) - - = 0, (A.ll) 

or r r z r 

^l + ^l + (Ag v + 2)W+Uh-h ff )=0, (A.12) 
dr r 2 

^ + ^ + (A ev + 2)X = 0. (A.13) 
dr r 

When decomposing the fields on a basis of spherical harmonics, these relations reduce 
to a system of ordinary differential equations with respect to r, which is solved by 
spectral methods in a similar way to the Poisson equation |25| . 
The numerical algorithm is then: 

(i) transform the source term of Eq. (|31|l to the Cartesian basis, 

(ii) solve the resulting six decoupled scalar Poisson equations for , 

(iii) transform back to the spherical basis and compute the potentials W and X, 

(iv) do an iteration on h, first solving the system i|A.ll|) - (|A.13|) with h, W and X as 
sources and then calculating the new value of h from the non-linear equation l|A.l|) . 

Since the system is overdetermined (four additional relations to satisfy), the 
integrability condition is that the source of the tensor Poisson equation 1|31JI be 
divergence-free. We do not impose this condition during the main iteration of the 
code, since this is not true for intermediate solutions of the metric and matter fields. 
We only check that this is satisfied, up to the accuracy of the code, at the end of 
the iteration, and we have found that this was true at the error level given by the 
virial identities (see Sec. 12. 4J) . The potentials W and X have been chosen among the 
six degrees of freedom of h %3 because none of their radial derivatives appear in the 
gauge conditions ljA.ll(l - l|A.13|) . hence, we do not calculate any radial derivative of 
these quantities to get the other components of h % i . Another reason for this choice is 
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that, when considering a more general case of dynamically evolving spacetime, these 
two potentials are asymptotically related to the two gravitational wave polarisation 
modes: P -> h + and h 9v -> h x in our asymptotically transverse-traceless gauge§. 
Note that in our case of stationary and axisymmetric spacetime, we have fj, = X = 0, 
which simplifies the resolution. 
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